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Abstract 

This paper describes the implementation of opti- 
mization techniques based on control theory for 
wing and wing-body design. In previous stud- 
ies [18, 19, 22] it was shown that control theory could 
be used to devise an effective optimization proce- 
dure for airfoils and wings in which the shape and 
the surrounding body-fitted mesh are both gener- 
ated analytically, and the control is the mapping 
function. Recently, the method has been imple- 
mented for both potential flows and flows governed 
by the Euler equations using an alternative formula- 
tion which employs numerically generated grids, so 
that it can more easily be extended to treat general 
configurations [34, 23]. Here results are presented 
both for the optimization of a swept wing using an 
analytic mapping, and for the optimization of wing 
and wing-body configurations using a general mesh. 

Nomenclature 

a spanwise scale factor 
A t Cartesian flux Jacobians 
A,B,C metric terms 

6 design variable 

6 bandwidth of full Jacobian matrix 
B generic boundary of flowfield domain 
Bw wing surface boundary of flowfield domain 
Bf far-field boundary of flowfield domain 
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Bs symmetry plane boundary of flowfield domain 
C\ transformed flux Jacobians 
C'p coefficient of drag 
Ci coefficient of lift 
C'p coefficient of pressure 
C* coefficient of pressure for sonic flow 
D flowfield domain 
E total energy 

ft components of Cartesian fluxes 

/, ft without pressure terms 

Ft components of transformed fluxes 

F x Ft without pressure terms 

T control law, physical location of boundary 

Go ^7,' 

Q gradient of design space 
Q projected Q into admissible space 
H total enthalpy 
/ cost function 

J Jacobian of generalized transformation 
K grid transformation matrix 
A. K at the profile surface 
h\j grid transformation matrix 
m number of design points 
M local Mach number 
A/oo Mach number at infinity 
n number of design variables 
h number of grid points 
fit components of a unit vector normal 
p pressure 
pd desired pressure 
V , Q , H metric terms 
q speed 

</oo speed at infinity 

R generic governing equation for flowfield 
normalized arc length 
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surface displacement in mapped plane 

reference area 

time 

exponents on basis functions 
Cartesian velocity components 
contravariant velocity components 
state vector of flowfield unknowns 
state vector scaled by J 
metric term 
mesh point positions 
Cartesian coordinates 
Cartesian coordinates 
normalized chord length 
wing surface coordinates 
Cartesian surface coordinates 
body fitted coordinates 
body fitted coordinates 

V + S 

modified boundary condition for V- 
angle of attack 
smoothing parameters 
variational 

Kronecker delta function 
ratio of specific heats 
distance along search direction 
smoothing parameters 
distance along search direction 
cost function weighting factors 
vector co-state variable 
density 

freestream density 


Introduction 

Since about 1960, there has been rapid progress in 
the field of Computational Fluid Dynamics (CFD). 
Especially in the last decade, with substantial im- 
provements in both computer performance and nu- 
merical methods, CFD has been extensively used 
together with experimental methods to aid in the 
aerodynamic design process. While much research 
continues in the CFD field, accurate and robust so- 
lutions are now routinely obtained over complete air- 
craft configurations for many flow conditions. Mod- 
em aircraft designers hope to benefit from this ca- 
pacity both to refine existing designs at transonic 
conditions and to develop new designs at supersonic 
conditions. These highly nonlinear flow regimes re- 
quire a design fidelity for which only CFD may pro- 
vide the answers within practical time constraints. 
Thus far however, CFD, like wind tunnel testing, ha & 
not had much success in direct aerodynamic shape 
design. Since the inception of CFD, researchers have 
sought not only accurate aerodynamic prediction 
methods for given configurations, but also design 
methods capable of creating new optimum config- 
urations. Yet, while flow analysis can now be car- 
ried out over quite complex configurations using the 
Navier-Stokes equations with a high degree of con- 


fidence, direct CFD based design is still limited to 
very simple two-dimensional and three-dimensional 
configurations, usually without the inclusion of vis- 
cous effects. The CFD-based aerodynamic design 
methods that do exist can be broken down into three 
basic categories: inverse surface methods, inverse 

field methods and numerical optimization methods. 

Inverse surface methods derive their name from 
the fact that they invert the goal of the flow analysis 
algorithm. Instead of obtaining the surface distribu- 
tion of an aerodynamic quantity, such as pressure, 
for a given shape, they calculate the shape for a 
given surface distribution of an aerodynamic quan- 
tity. An alternative way to obtain desirable aero- 
dynamic shapes is through the use of field-based in- 
verse design methods. These methods differ from 
surface specification methods in that they obtain 
designs based upon objectives, or constraints, im- 
posed not only upon the configuration surface but 
everywhere in the flow field. For transonic flows a 
field-based objective can be to limit shock strength 
or create shock-free designs. Most of these methods 
are based on potential flow techniques, and few of 
them have been extended to three-dimensions. 

The common trait of all inverse methods is their 
computational efficiency. Typically, transonic in- 
verse methods require the equivalent of 2-10 com- 
plete flow solutions in order to render a complete 
design. Since obtaining a few solutions for simple 
two-dimensional and three-dimensional designs can 
be done in at most a few hours on modern comput- 
ers systems, the computational cost of most inverse 
methods is considered to be minimal. Unfortunately, 
they suffer from many limitations and difficulties. 
Their most glaring limitation is that the objective is 
built directly into the design process and thus cannot 
be changed to an arbitrary or more appropriate ob- 
jective function. The user must therefore be highly 
experienced in order to be able to prescribe surface 
distributions or choose initial geometries which lead 
to the desired aerodynamic properties. In addition, 
surface inverse methods have a tendency to fail be- 
cause the target surface distribution is not necessar- 
ily attainable. In general it must satisfy constraints 
to permit the existence of the desired solution. On 
the other hand, field inverse methods typically only 
allow for the design of a single shock-free design 
point, and have no means of properly addressing 
off-design conditions. Furthermore, it is difficult to 
formulate inverse methods that can satisfy desired 
aerodynamic and geometric constraints. In essence, 
inverse methods require designers to have an a pri- 
on knowledge of an optimum pressure distribution 
that satisfies the geometric and aerodynamic con- 
straints. This limited design capability and difficult 
implementation has, to date, constrained the appli- 
cability of inverse methods. 

An alternative approach, which avoids some of the 
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difficulties of inverse methods, but only at the price 
of heavy computational expense, is to use numerical 
optimization methods. The essence of these meth- 
ods is very simple: a numerical optimization pro- 
cedure is coupled directly to an existing CFD anal- 
ysis algorithm. The numerical optimization proce- 
dure attempts to extremize a chosen aerodynamic 
measure of merit which is evaluated by the chosen 
CFD code. The configuration is systematically mod- 
ified through user specified design variables. These 
design variables must be chosen in such a way as 
to permit the shape of the configuration to change 
in a manner that allows the design objective to be 
improved. Most of these optimization procedures 
require gradient information in addition to evalua- 
tions of the objective function. Here, the gradient 
refers to changes in the objective function with re- 
spect to changes in the design variables. The sim- 
plest method of obtaining gradient information is by 
“brute-force” finite differences. In this technique, 
the gradient components are estimated by indepen- 
dently perturbing each design variable with a finite 
step, calculating the corresponding value of the ob- 
jective function using CFD analysis, and forming 
the ratio of the differences. The gradient is used 
by the numerical optimization algorithm to calcu- 
late a search direction using steepest descent, con- 
jugate gradient, or quasi-Newton techniques. The 
optimization algorithm then proceeds by estimating 
the minimum or maximum of the aerodynamic ob- 
jective function along the search direction using re- 
peated CFD flow analyses. The entire process is re- 
peated until the gradient approaches zero or further 
improvement in the aerodynamic objective function 
is impossible. 

The use of numerical optimization for transonic 
aerodynamic shape design was pioneered by Hicks, 
Murman and Vanderplaats [13]. They applied the 
method to two-dimensional profile design subject to 
the potential flow equation. The method was quickly 
extended to wing design by Hicks and Henne [11, 12] . 
Recently, in the work of Reuther, Cliff, Hicks and 
Van Dam, the method has proven to be successful for 
the design of supersonic wing/body transport config- 
urations through its extension to three-dimensional 
flows governed by the Euler equations [33]. In all 
of these cases brute-force finite difference methods 
were used to obtain the required gradient informa- 
tion. 

These methods are very versatile, allowing any 
reasonable aerodynamic quantity to be used as the 
objective function. They can be used to mimic an in- 
verse method by minimizing the difference between 
target and actual pressure distributions, or may in- 
stead be used to maximize other aerodynamic quan- 
tities of merit such as L/D. Geometric constraints 
can be readily enforced by a proper choice of design 
variables. Aerodynamic constraints can be treated 


either by adding weighted terms to the objective 
function or by the use of a constrained optimization 
algorithm. Unfortunately, these brute-force numer- 
ical optimization methods, unlike the inverse meth- 
ods, are computationally expensive because of the 
large number of flow solutions needed to determine 
the gradient information for a useful number of de- 
sign variables. For three-dimensional configurations, 
hundreds or even thousands of design variables may 
be necessary. This implies that tens of thousands 
of flow analyses would be required for a complete 
design. 

Formulation of the design 
problem as a control problem 

Clearly, alternative methods must be developed 
which have the flexibility and power of current nu- 
merical optimization codes but do not require such 
large computational resources. These new methods 
must avoid the limitations and difficulties of tradi- 
tional inverse methods while approaching their in- 
herent computational efficiency. 

One means of attaining such a design method is 
by treating the design problem within the mathe- 
matical theory for control of systems governed by 
partial differential equations [29]. Suppose that the 
boundary is defined by a function T( b ), where b is 
the position vector of the design variables, and the 
desired objective is measured by a cost function /. 
This may, for example, measure the deviation from 
a desired surface pressure distribution, but it can 
also represent other measures of performance such 
as the drag or the lift/drag ratio. Suppose that a 
variation ST in the control produces a variation SI 
in the cost. SI can be expressed to first order as an 
inner product 

SI = (£, ST ) , 

where the gradient, £/, of the cost function with re- 
spect, to the control, is independent of the particu- 
lar variation ST . Following control theory Q can be 
determined by solving an adjoint equation. If one 
makes a shape change 

ST = —AC/, 

where A is sufficiently small and positive, then 

SI = - A(S,£)<0 (1) 

assuring a reduction in /. The method can be accel- 
erated by choosing ST not simply as a multiple of 
the gradient (steepest descent) but instead as a more 
sophisticated search direction provided by numerical 
optimization. 

For flow about an airfoil or wing, the aerodynamic 
properties which define the cost function are func- 
tions of the flow-field variables ( te), the physical lo- 
cations of the mesh points within the volume (,V), 
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( 5 ) 


and the physical location of the boundary (T). Then 


the first term is eliminated, and we find that 


7 = I (w, X , T) , 

and a change in T results in a change 
df T d! r dI T 

4 ' = isr'" , ’ + a x tx + W 6r ' (2) 

in the cost function. As pointed out by Baysal and 
Eleshaky [4] each term in (2), except for Sw , can be 
easily obtained. ^ and jj can be obtained di- 
rectly without a fiowfield evaluation since they are 
partial derivatives. ST is simply the surface modifi- 
cation and SX can be determined by either working 
out the exact analytical values from a mapping, or 
by successive grid generation for each design vari- 
able, so long as this cost is significantly less then 
the cost of the flow solution. For solutions requiring 
a large number of mesh points where grid genera- 
tion becomes expensive, an alternative method for 
calculating SX can be formulated using grid pertur- 
bation. Brute- force methods evaluate the gradient 
by making a small change in each design variable 
separately, and then recalculate both the grid and 
flow-field variables. This requires a number of addi- 
tional flow calculations equal to the number of de- 
sign variables. Using control theory, the governing 
equations of the fiowfield are introduced as a con- 
straint. in such a way that the final expression for the 
gradient does not require multiple flow solutions. In 
order to achieve this result, Sw must be eliminated 
from (2). The governing equation R expresses the 
dependence of w, X and T within the fiowfield do- 
main D, 

R(w, X , T) = 0, 


Thus is determined from the equation 


SR 


OR 

dw 


Sw + 


dR 

dx 


sx + 


dR 

dr 


ST- 0. (3) 


Next, introducing a Lagrange Multiplier t p, we have 
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Choosing i /’ to satisfy the adjoint equation 


dR 


dw 


, 81 


(4) 


61 = G t ST 

where 



(6) 


The advantage is that (5) is independent of Sw, with 
the result that the gradient of 7 with respect to 
an arbitrary number of design variables can be de- 
termined without the need for additional flow-field 
evaluations. The main cost is in solving the ad- 
joint equation (4). In general, the adjoint problem is 
about as complex as a flow solution. However, if the 
number of design variables is large, the cost differen- 
tial between one adjoint solution and the large num- 
ber of fiowfield evaluations required to determine the 
gradient by finite differences becomes compelling. 

Once the gradient is calculated by (6) a modifica- 
tion following (1) can then be made. After making 
such a modification, the gradient can be recalcu- 
lated and the process repeated to follow a path of 
steepest descent until a minimum is reached. In or- 
der to avoid violating geometric constraints, such 
as a minimum acceptable wing thickness, the gra- 
dient may be projected into the allowable subspace 
within which the constraints are satisfied. In this 
way procedures can be devised which must necessar 
ily converge at least to a local minimum. The effi- 
ciency may be improved by performing line searches 
to find the minimum in a search direction defined 
by the negative gradient, and also by the use of 
more sophisticated descent methods such as conju- 
gate gradient or quasi-Newton algorithms. There 
is the possibility of more than one local minimum, 
but in any case the method will lead to an improve- 
ment over the original design. Furthermore, unlike 
the traditional inverse algorithms, any measure of 
performance can be used as the cost function. 

In this research, the adjoint approach is used 
to permit a dramatic reduction in the computa- 
tional cost of each design solution, since the gradi- 
ent cost will be reduced to the cost of approximately 
two flow evaluations (provided the adjoint equations 
are about as computationally expensive as the flow 
equations) instead of the traditional n + 1 evalua- 
tions, where n is the number of design variables. 
The key point is that the cost of the optimization 
method is no longer proportional to the number of 
design variables, which has been the limiting factor 
in brute-force aerodynamic optimization methods. 

Another significant advantage of this method is 
its applicability to multi-point design problems. Be- 
cause of its efficient use of computer resources, two 
or three design points can be included in the opti- 
mization procedure by solving separate adjoint prob- 
lems for each design point. The resulting procedure 
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yields a method which is capable of calculating the 
combined gradient from all design points from order 
2 m equivalent flow calculations, with m being the 
number of design points. Even in the case where 
there is a large number of design points a significant 
benefit is still realized when compared with brute- 
force calculations where the number of solutions re- 
quired is mn. The use of this approach will allow the 
performance benefits at various design points to be 
considered together, thus obtaining a more optimal 
overall design. 

A variety of alternative formulations of the design 
problem can then be treated systematically within 
the framework of the mathematical theory for con- 
trol of systems governed by partial differential equa- 
tions [29]. The use of this method for aerodynamic 
design was first introduced by Jameson [18] who ex- 
amined its application to transonic flow governed by 
both the potential flow and Euler equations [19, 22]. 
In his work, control theory was applied directly to 
the partial differential equations governing the flow 
solution. Thus the adjoint equations were formed as 
a system of differential equations. These adjoint dif- 
ferential equations were then discretized and solved 
in the same manner as the flow equations to obtain 
the necessary gradient information. This approach, 
often termed the continuous sensitivity analysis, was 
used by Jameson in conjunction with analytic grid 
mappings, to formulate directly at the differential 
level, the necessary systems of equations defined by 
equation ( 6 ). 

Steps (2 - 6 ) may alternatively be applied to 
the discrete equations which approximate the gov- 
erning differential equations. This approach (dis- 
crete sensitivity analysis) is now gaining favor in the 
work of Korovi, Newman, Taylor, Hou and Jones 
[30, 26, 15, 24] and also Baysal, Eleshaky and Bur- 
green [3, 4, 5, 7, 2, 8 , 9, 6 ]. The continuous and 
the discrete formulations methods can be very sim- 
ilar, depending upon the discretization of (4). The 
method used by Jameson has an advantage in that 
the discretization and iteration scheme used to solve 
the flowfield system can also be applied directly to 
the adjoint system. Therefore, the robust iteration 
algorithms and convergence acceleration techniques 
that have been matured for CFD algorithms can be 
directly ported for the solution of the adjoint sys- 
tem. Discrete sensitivity methods for equations (2 
- 6 ) often resort to matrix elimination methods to 
solve (4). While direct techniques to solve these 
large sparse systems can be robust and reliable they 
suffer when the number of points becomes large be- 
cause the operation count grows as 0(nb 2 ) and the 
storage goes as 0 (n 6 ), where n is the number of 
unknowns and b is the bandwidth. Therefore, in or- 
der to solve these large systems, alternatives such 
as sophisticated matrix decomposition [28] or in- 
cremental iterative methods [25] must be employed. 


However, these alternatives, whether applied to the 
flow solution system or, as here, to the adjoint sys- 
tem, have not proven to be as efficient as meth- 
ods developed for CFD. CFD approaches altogether 
avoid the direct matrix elimination process by re- 
laxing the system either point-by-point or line-by- 
line in a pseudo-time procedure. When the num- 
ber of mesh points becomes large, especially in the 
case of three-dimensional problems, the 0{h) oper- 
ational counts and the O(h) storage of explicit iter- 
ation schemes used in CFD can significantly reduce 
the time and memory required to solve the adjoint 
system. Methods that use matrix decomposition or 
the Generalized Minimum Residual (GMRES) iter- 
ative approach to solve the discrete sensitivity prob- 
lem have shown modest improvements over standard 
Gaussian elimination [28] and [37, 14]. Nevertheless, 
these methods are currently still not competitive in 
terms of both time and memory with mature CFD- 
like algorithms. The efficiency of discrete sensitivity 
methods may be improved by the introduction of in- 
cremental iterative strategies for solving the matrix 
elimination problem. Present work by Taylor et al. 
[25] to develop such methods shows great promise. 
Currently the use of continuous sensitivity analy- 
sis eliminates the need for developing incremental 
strategies since the existing (and presumably ma- 
ture) flow solution algorithm can be used directly 
for the adjoint solution. 

Other applications of control theory in aerody- 
namics have been explored by Pironneau for opti- 
mum shape design of systems governed by elliptic 
equations [31]. More recently Ta’asan, Kuruvila, 
and Salas implemented a one shot approach in which 
the constraint represented by the flow equations is 
only required to be satisfied by the final converged 
solution [36, 27]. Pironneau has also studied the use 
of control theory for optimum shape design of sys- 
tems governed by the Navier-Stokes equations [32]. 


Three-Dimensional Design 
Using the Euler Equations 

The application of control theory to aerodynamic 
design problems is illustrated by treating the case 
of three-dimensional wing design, using the inviscid 
Euler equations as the mathematical model for com- 
pressible flow. In this case it proves convenient to 
denote the Cartesian coordinates and velocity com- 
ponents by a*i, x 2 , £3 and u i, 1 / 2 , u 3 , and to use the 
convention that summation over i = 1 to 3 is im- 
plied by a repeated index i. The three-dimensional 
Euler equations may be written as 


5 



where 


\ 

p 


f \ 

pu, 

pu\ 


pu,u x + p6 ti 

pu 7 

) , /. = < 

pu,u 2 + pb, 2 

pU 3 


pu,u 3 + pbi 3 

PE J 


„ puiH 


and b tJ is the Kronecker delta function. Also, 

P= (7 - Opjtf- ^ («<)| , (9) 


and 

pH = pE + p ( 10 ) 

where 7 is the ratio of the specific heats. Consider 
a transformation to coordinates £j , £3 where 


Ki, 



J = det (H) , h'r. 1 



Introduce contravariant velocity components as 


As the first example, suppose that it is desired 
to control the surface pressure by varying the wing 
shape. It is convenient to retain a fixed computa- 
tional domain. Variations in the shape then result 
in corresponding variations in the mapping deriva- 
tives defined by K. Introduce the cost function 

1 = \ J j B 

where pd is the desired pressure on the wing. The 
design problem is now treated as a control prob- 
lem where the control function is the wing shape, 
which is to be chosen to minimize / subject to the 
constraints defined by the flow equations ( 7 - 12 ). A 
variation in the shape will cause a variation bp in 
the pressure and consequently a variation in the cost 
function 


bl = 



{p~Pd)bp d^d^ 3 . 


( 14 ) 
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The Euler equations can now be written as 


with 


W 


dW dF t 

~dT + ^ 
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plhu 2 + M;P > 

PU 3 


pUiU 3 + §^p 

pE j 


pUiH 
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(12) 


Since p depends on w through the equation of 
state ( 9 - 10 ), the variation bp can be determined 
from the variation bw. Define the Jacobian matrices 

A ‘ = %- <■ os) 

Then the equation for bw in the steady state be- 
comes 

where 

6Ft =Ci6w+6^ji^J fj. 

Now, multiplying by a vector co-state variable V’ and 
integrating over the domain 



Assume now that designs are limited to wing and 
wing-body configurations using a C-H topology 
mesh. If a body is present it is assumed to exist 
as part of the face containing the symmetry plane. 
Thus the new computational coordinate system con- 
forms to the wing in such a way that the wing surface 
Bw is represented by £2 = 0 and the body and sym- 
metry plane B$ conform such that £3 = 0 . Then 
the flow is determined els the steady state solution 
of equation (11) subject to the flow tangency condi- 
tions 


[l 2 = 0 on Bw 

U 3 = 0 on 0 5 . (13) 

At the far field boundary Bp, freestream conditions 
are specified for incoming waves, as in the two- 
dimensional case, while outgoing waves are deter- 
mined by the solution. 


and if V’ is differentiable this may be integrated by 
parts to give 

where h, are components of a unit vector normal to 
the boundary. Thus the variation in the cost func- 
tion may now be written 
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On the wing surface Bw , »i = n 3 = 0 and on the 
body and symmetry plane B. s , fi , = fi 2 = 0 . It 



follows from equation (13) that 
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Suppose now that F is the steady state solution 
of the adjoint equation 




(18) 


At the outer boundary incoming characteristics for 
t /’ correspond to outgoing characteristics for 6w. 
Consequently, as in the two-dimensional case, one 
can choose boundary conditions for V- such that 


iiiil? CiSw = 0. 

Then if the coordinate transformation is such that 
£ (J/\ _1 ) is negligible in the far field, the only re- 
maining boundary terms are 

-// tl> T 6Fi dZid(,3, 

J JBw 

- JJ rJ' T 6F 3 d^d^ 2 . 

Thus by letting t/' satisfy the boundary conditions, 

J (v' 2 ^- + 03“^- + = {p ~ Pd) on Bw, 

J 4* ^4^- + Vm 4^-) = 0 on 0 9 ) 

\ aii ox 2 oxz/ 

we find finally that 

{ v2< ( j £f) 4 ( J S) + ( J S) 
{* 2> ( j ^) + *•** ( j ^t) } 


^ ^ B , 

JL 


s? " u 
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la x,y-Plane. 

lb: r;* Plane. 

Figure 1: Sheared Parabolic Mapping. 


Wing Design Using Analytic 
Mesh Generation 

A convenient way to treat a wing is to introduce 
sheared parabolic coordinates as shown in figure 1 
through the transformation 

* = *o(C) + lo(C) |^ 2 - (»/+ S(£,C)) 2 } 

y = vo (O + aKKto + ^K.O) 

- = C- 

Here j = n, y = z = x 3 are the Cartesian 
coordinates, and £ and correspond to parabolic 

coordinates generated by the mapping 

1 2 

X + iy - *0 + iyo + -a (<) {(, + if)} 

at a fixed span station £ and ij = r; + .S’, xo (C) and 
y 0 (C) are the coordinates of a singular line which is 
swept to lie just inside the leading edge of a swept 
wing, while a (£) is a scale factor to allow for span- 
wise chord variations. The surface r/ = 0 is a shal- 
low bump corresponding to the wing surface, with a 
height S (£, C) determined by the equation 


t + = ^(xbvv + iVB w ), 

where xb w (-) and ys w ( 2 ) are coordinates of points 
lying on the wing surface. We now treat S (£, C) as 
the control. 

In this case the transformation matrix 4?* be- 
comes 
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Then under a modification bS 

bx^ — — a (SSSf + i)6S{) 

bx v — — aSS 

hz = <i(6S + £6St) 

by v ~ 0. 


Thus 

and 
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0 

sj 


where 


W = 6y e .4 - bx ( B - a c --65 - 6J5< - JbS 2 

Inserting these formulas in equation (20) we find 
that the volume integral in 61 is 


/// 


d\ 

~^—bSf 2 dt) dC 
dij )T 


bl= I I (r(SX)bS-Q(S<C)Mt-K(tX)6S<)dtdc 

' Bu 


- J J I + 6x ^ 2 + di) ^ 

+ IIJ 3 K d n 

where $ and bS are independent of //. Therefore, 
integrating over 7/, the variation in the cost function 
can be reduced to a surface integral of the form 

IL 

Here 

T = a(i’ 2 + S( V’3 + CM p 

~ J {Zf 1 + vf 2 + (M + nB) h) dr] 

f dt T 

- J -jj-{h+S(f2+Ch)d V 

-Jw jd ” 

Q = a ((v’2 + vV’3 ) V 

dt T 


J 


dr, 


itf 1 + '7/2 + ((.A + i)B) h) di) 


7v = Jli'4p 


+ J ^J(’Adi), 


(21) 


C = 2af)S ( - A - BS ( + 

a 

Also the shape change will be confined to a boundary 
region of the £ plane, so we can integrate by parts 
to obtain 





Thus to reduce I we can choose 


65 = -A 


< dQ 

\ <K 



where A is sufficiently small and non-negative. 

In order to impose a thickness constraint we 
can define a baseline surface S 0 (£>C) below which 
5(£,C) is not allowed to fall. Now if we take 
A = A(£,C) as a non-negative function such that 


S(€,C) + *SK,C) >So(€,0- 


Then the constraint is satisfied, while 




2 

d£ dC < 0 


or 


c =' ,+ f 



( 22 ) 


To make sure that no discontinuities are introduced 
in the modified shape, smoothing can be introduced 
in the update process by setting 


0i6Z-^-02^zSZ = 

(>(. 


= -(' 


-f + 



(23) 


Then, integrating by parts in the <£ direction gives 


bl 





ihbZ- 


d ^ 



dt) 


Wing-Body Design Using 
a General Mesh 


In order to treat a more general mesh we revert 
to equations (14-20). The difficulty in using these 
equations is that the variation of the metric terras 
in the equations needs to be obtained in order to 
construct bl in equation (20). One way to accom- 
plish this is to use finite differences to calculate the 
necessary information. While this approach would 
avoid the use of multiple flow solutions to determine 
the gradient, it would unfortunately still require the 
mesh generator to be used repeatedly. The number 
of mesh solutions required would be proportional to 
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the number of design variables. This may be accept- 
able, since the flow solution process is typically much 
more computationally expensive than grid genera- 
tion. Such a method should then ensure a significant 
savings over using finite differences for both the grid 
generation and flow solution processes. However, 
for three-dimensional designs where both the num- 
ber of design variables and the computational cost 
of grid generation can be high, this method is ex- 
cessively expensive. Further, for complicated three- 
dimensional configurations, for which it is still not 
practical to integrate fully automatic grid genera- 
tion into the solution process, the method will not 
be feasible. 

This motivates the need to find a method which 
by-passes these difficulties. In order to remove 
the cost of the successive grid generation from the 
gradient calculation, a successive grid perturbation 
method is therefore used. In this approach, which 
was also used by Burgreen and Baysal [7], an initial 
structured curvilinear body-fitted grid over the ini- 
tial configuration is created by any grid generation 
process before optimization. Then the geometry as 
well as the grid become inputs to the optimization 
process. New grids, which conform to the surface 
as it is modified, can then be generated by shifting 
the grid points along each grid index line projecting 
from the surface by an amount which is attenuated 
as the arc length from the surface increases. If the 
outer boundary of the grid domain is held constant 
the modification to the grid has the form 


neu, = old 


+ S otd (x™ w - x°[ d ) , 


where represents the volume grid points, x 3x rep- 
resents the surface grid points and 5 represents the 
arclength along the radial mesh line measured from 
the outer domain, normalized so that S = 1 at the 
inner surface. The required variations in the met- 
ric terms can then be obtained in terms of surface 
perturbations since, 


6xi = S oU 6x 3 '. 


and 

f Kj 

We introduce, 


(24) 


G„ = J K ~ 1 



' VrjZ C - y< Z V 

X £ Zrj 

- XcVr, 

= 

y< z t - vt*< 

x<z< - x c z* 

- x *y< 


. ViZr) - y v z i 

Xr \ z i ~ \ 

Vv ~ . 


Now it is convenient to rewrite equation (20) after 
integration by parts as 





4> T 6G 2j fj } pd^idfr 



\1> T SGzjfj } p d£\ (ff 2 , 


( 25 ) 


where fj represent the flux components fj with the 
pressure terms dropped. From the definition of G tJ 
we have, for example, 


f>G n = S(yrj) z< + 6 (z<;)y n - S (j/<) z v - S (z„) y<; 

= $° ld [* (yr), ) ) Vn - b {y<. ) z n - 6 (zr). ) y<] . 

( 26 ) 


Substituting these expressions into equation (25) al- 
lows us to integrate along the index direction pro- 
jecting from the configuration surface without any 
dependence on particular design variables, since the 
metric variations are fully determined by the surface 
perturbations. Thus, the expression for the varia- 
tion in the cost function can be reduced to surface 
integrals only. 

While this type of grid perturbation method does 
not guarantee that grid lines will not eventually 
cross if the perturbations are large, this point is 
irrelevant for gradient calculations since only ana- 
lytic grid derivative information is needed. However, 
since we employ a numerical optimization algorithm 
with line searches along a descent direction, true re- 
gridding is also necessary. For these line search cal- 
culations the grid perturbation algorithm is used so 
long as negative cell volumes are not created. If sin- 
gularities begin to develop in the grid, the original 
grid generator can be used to create a new grid and 
the process restarted. In this work a modified ver- 
sion of WBGRID is used for automatic generation 
of the base grid, and subsequent reinitializations of 
the grid if needed. 

After substitution of (26), the resulting expres- 
sion for SI is reduced to surface integrals in which 
the remaining unknowns are the grid metrics. These 
surface grid metrics can be easily determined for 
any modification in the surface by direct evalua- 
tion. This suggests choosing a set of design variables 
which smoothly modifies the original shape, say 6,. 
The gradient can then be defined with respect to 
these design variables as 

G{b ' ) = ji' (27) 

where SI is calculated by (25) with each term b t 
being independently perturbed by a finite step. 

Therefore, to construct Q , a basis space of inde- 
pendent perturbation functions 6 t , i — l,2,...,n 
(n = number of design variables) should be cho- 
sen that allows for the needed freedom of the design 
space. In this work design variables have been cho- 
sen with the following chord wise form, suggested by 
Hicks and Henne [11, 12]: 


b(x) = sin ( 7r x i o i ) 


<2 
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b ( x ) = x ix ( 1 — x)e ta# , 

where t\ and 1 2 control the center and width of the 
perturbation and x is the normalized chord length. 
When distributed over the entire chord on both up- 
per and lower surfaces these analytic perturbation 
functions admit a large design space for each wing 
defining station. Then, by choosing a number of 
defining stations which reflects a practical design 
and linearly lofting the spanwise variations, the en- 
tire wing may be designed. An additional twist de- 
sign variable is also included at each station to allow 
for washout that may be needed. These design vari- 
ables have the advantage of being space-based func- 
tions, as opposed to frequency-based functions, and 
thus they allow for local control of the design. They 
can be chosen such that symmetry, thickness, or vol- 
ume can be implicitly constrained. Further, partic- 
ular choices of these variables will concentrate the 
design effort in regions where refinement is needed, 
while leaving the rest of the wing virtually undis- 
turbed. The disadvantage of these functions is that 
they do not necessarily form a complete basis space 
when their number is increased, nor are they orthog- 
onal. Thus, they do not guarantee that a solution, 
for example, of the inverse problem for a realizable 
target pressure distribution will necessarily be at- 
tained. Here, they are employed for their ease of 
use and ability to produce a wide variation of shapes 
with a limited number of design variables. Note that 
th is limited set of design variables does not allow 
for planform modifications. However, here the in- 
terest is focused on Aerodynamic optimization and 
not on true Multi-Disciplinary optimization. In or- 
der to accomplish realistic planform design a true 
multi-point multi-disciplinary design method incor- 
porating all aspects of aircraft direct operating cost 
would have to be used, because of trade-offs such as 
that between wing span and weight. Here the prob- 
lem focuses on squeezing the most performance from 
a given planform. 

Implementation of the Euler based 
design methods 

Both of the design methods have been successfully 
implemented. The two techniques share many com- 
mon features such as the flow and adjoint solution 
algorithms. The procedures can be summarized as 
follows. 

1. Solve the flow equations (7—12) for p , u t , p, £, 
//, and (/, . 

2. Smooth the cost function if necessary by (23). 

3. Solve the adjoint equations (18) for V ; subject 
to the boundary conditions (19). 


4. Either calculate V , Q and 72., by (21), from the 
variation in the control £(£,(), or evaluate the 
surface independent terms in equation (26). 

5. Evaluate Q by equation (22) or (27). 

6. Then either correct the mapping S (£,£) or up- 
date the design variables fr t based on the direc- 
tion from steepest descent 

6S {CO = -A(? or Sb t = -XQ 
or as an alternative a quasi-Newton method. 

7. Return to 1. 

In practice these methods resemble those used by 
Hicks, Reuther et al. [35, 13, 33] with control the- 
ory replacing the brute-force, finite difference gradi- 
ent calculation. Unlike the earlier procedures, the 
current methods’ computational costs do not hinge 
upon the number of design variables, which in these 
cases is either the number of surface mesh points 
used to represent S (£,£), or the number of perturba- 
tion functions b t . Thus, with the three-dimensional 
implementation in hand, nonlinear aerodynamic de- 
sign of complete aircraft can be brought into the 
realm of computational feasibility. The method also 
has the advantage of being quite general in that ar- 
bitrary choices for both the design variables and op- 
timization technique are admitted. 

The practical implementation of the design meth- 
ods rely heavily upon fast and accurate solvers for 
both the state (w) and co-state (V\) fields. Fur- 
ther, to improve the speed and realizability of the 
methods, a robust choice of the optimization algo- 
rithm must be made. In the present case the sec- 
ond author’s FL087 computer program has been 
used as the basis of the design method FL087 
solves the three-dimensional Euler equations with 
a cell-centered finite volume scheme, and uses resid- 
ual averaging and multigrid acceleration to obtain 
very rapid steady state solutions, usually in 25 to 
50 multigrid cycles [16, 17]. Upwind biasing is used 
to produce nonoscillatory solutions, and assure the 
clean capture of shock waves. This is introduced 
through the addition of carefully controlled numeri- 
cal diffusion terms, with a magnitude of order Ax 1 2 3 
in smooth parts of the flow'. The adjoint equations 
are treated in the same way as the flow equations. 
The fluxes are first estimated by central differences, 
and then modified by downwind biasing through nu- 
merical diffusive terms. These are supplied by the 
same subroutines as are used for the flow equations. 

In the implementation with analytic mapping, a 
gradient procedure is used as the optimization pro- 
cess. However, to preserve the smoothness of the 
profile the gradient is smoothed at each step. Thus 
the change in the shape function S (£,0 is defined 
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by solving 


6S -TZl i T£ 6 S =-M' 


where (3 is a smoothing parameter. Then, to first 
order, the variation in the cost is 


61 


< 



Thus an improvement is still asured when smoothing 
is used. For the implementation on arbitrary meshes 
a quasi-Newton optimization method is employed. 
For this purpose the QNMDIF program developed 
by Gill, Murray and Wright [10] is used. 

The option to minimize the pressure drag coeffi- 
cient is realized in both methods by redefining the 
cost function as 

/ = C D = — f [ pG 2 id£id£ 3 

2 Poo<7oo ^ref ^ JBw 

+ 1 f / pC*3i d £ i d^2 

2Poo9oo^ ref J JBs 

where S re f is the reference area. To prevent the pro- 
cedure from trying to reduce drag by reducing the 
profile to a non-lifting flat plate a target pressure dis- 
tribution can be retained in the cost function, which 
becomes 

/ = ^fii f f ( p-Pd)bp d£\d (, 3 + 

*■ J JBw 

where and are weighting parameters. Similar 
constructions are employed for other cost functions 
such as (1/-^). 


Numerical tests of the 
Three-Dimensional Method 
Using an Analytic Mapping 

The analytic grid generation method has been tested 
for the optimization of a swept wing [22, 21]. Two 
examples are presented here. In each the planform 
was fixed while the wing sections were free to be 
changed arbitrarily by the design method. 

In the first example the wing has a unit-semi- 
span, with 36 degrees leading edge sweep. It has 
a compound trapezoidal planform, with straight ta- 
per from a root chord of 0.38 to a chord of 0.26 at 
the 30 percent span station, and straight taper from 
there to a chord of 0.12 at the tip, with an aspect 


ratio of 8.7. The initial wing sections were based 
on the Korn airfoil, which was designed for shock 
free flow at Mach 0.75 with a lift coefficient of 0.63, 
and has a thickness to chord ratio of 11.5 percent 
[1]. The thickness to chord ratio was increased by 
a factor of 1.2 at the root and decreased by a ratio 
of 0.8 at the tip, with a linear variation across the 
span. The inboard sections were rotated upwards to 
give 3.5 degrees twist across the span. 

The two-dimensional pressure distribution of the 
Korn airfoil at its design point was introduced as 
a target pressure distribution uniformly across the 
span. This target is presumably not realizable since 
it would imply a lifting wing with zero vortex drag, 
but it serves to favor the establishment of a rela- 
tively benign pressure distribution. The total invis- 
cid drag coefficient, due to the combination of vortex 
and shock wave drag, was also included in the cost 
function. Calculations were performed with the lift 
coefficient forced to approach a fixed value by ad- 
justing the angle of attack every fifth iteration of 
the flow solution. A grid with 192x32x48 = 294,912 
points was used, and the wing shape was determined 
by 133 sections each with 128 mesh points for a to- 
tal of 4,224 design variables. It was found that the 
computational costs can be reduced by using only 
15 multigrid cycles in each flow solution, and in each 
adjoint solution. Although this is not enough for full 
convergence, it proves sufficient to provide a shape 
modification which leads to an improvement. Fig- 
ures 2, 3, and 4 shows the result of a calculation at 
Mach number of 0.82, with the lift coefficient forced 
to approach a value of 0.5. The plots show the initial 
wing geometry and pressure distribution, and the 
modified geometry and pressure distribution after 8 
design cycles. The total inviscid drag was reduced 
from 0.0185 to 0.0118. The initial design exhibits a 
very strong shock wave in the inboard region. It can 
be seen that this is completely eliminated, leaving a 
very weak shock wave in the outboard region. 

To verify the solution, the final geometry, after 8 
design cycles, was analyzed with another method, 
using the computer program FL067. This program 
uses a cell- vertex formulation, and has recently been 
modified to incorporate a local extremum diminish- 
ing (LED) algorithm with a very low level of numer- 
ical diffusion [20]. When run to full convergence it 
was found that the initial wing has a drag coeffi- 
cient of 0.0171 at Mach 0.82 and a lift coefficient of 
0.496, with a corresponding lift to drag ratio of 29. 
The result is illustrated in Figure 5. The final wing, 
shown in Figure 6, has a drag coefficient of 0.0107 
with a lift coefficient of 0.497, giving a lift to drag 
ratio of 47 at the same Mach number. A calculation 
at Mach 0.50 shows a drag coefficient of 0.0100 for 
a lift coefficient of 0.500. Since in this case the flow 
is entirely subsonic, this provides an estimate of the 
vortex drag for this planform and lift distribution. 
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Thus the design method has reduced the shock wave 
drag coefficient to about 0.0007. 

The second example is a design at a higher Mach 
number of 0.85, with a target lift coefficient of 0.55. 
This is a more severe test of the method, and a 
higher sweep back angle of 38 degrees at the leading 
edge was used. The wing has a modified trapezoidal 
planform, with straight taper from a root chord of 
0.38, and a curved trailing edge in the inboard re- 
gion blending into straight taper outboard of the 30 
percent span station to a tip chord of 0.10, with an 
aspect ration of 9.0. The initial wing sections were 
based on a section specially designed by the second 
author’s two dimensional design method [18] to give 
shock free flow at Mach 0.78 with a lift coefficient 
of 0.6. This section, which has a thickness to chord 
ratio of 9. 5 percent, was used at the tip. Similar sec- 
tions with an increased thickness were used inboard. 
The variation of thickness was nonlinear with a more 
rapid increase near the root, where the thickness to 
chord ratio of the basic section was multiplied by 
a factor of 1.47. The inboard sections were rotated 
upwards to give the initial wing 3 degrees of twist 
from root to tip. The two dimensional pressure dis- 
tribution of the basic wing section at its design point 
again was introduced as a target pressure distribu- 
tion uniformly across the span. Figures 7 and 8 
show the result of the optimization. The total in- 
viscid drag coefficient was reduced from 0.0243 to 
0.0144. The results of this optimization were also 
verified by calculations with FL067, using a high 
resolution LED algorithm. Figures 9 and 10 show 
that when the solutions were fully converged the 
drag coefficient was reduced from 0.0236 to 0.0119, 
with an improvement in lift to drag ratio from 23 to 
46. The result is illustrated in figure 10. A subsonic 
calculation at Mach 0.50 shows a drag coefficient of 
0.0109 for the same lift coefficient of 0.55. Thus in 
this case the shock wave drag coefficient is about 
0.0010. 

Numerical tests of the 
Three-Dimensional Method 
for Wing and Wing-Body 
Configurations Using a 
General Mesh 

The first design case involving the arbitrary mesh 
implementation is an inverse design to obtain the 
ONERA-M6 wing at M = 0.84 and a = 3.06°. The 
fixed alpha design starts from the ONERA-M6 plan- 
form but has NACA 0012 airfoil sections instead of 
the original sections. The goal is thus to recover the 
original ONERA-M6 sections by prescribing its ac- 
tual pressure distribution at the desired conditions 
as the target. Design variables are specified at 6 
defining stations which are then lofted in the span- 


wise direction. Twenty-five design variables are used 
at each defining station, for a total of 150. They are 
specified such that only thickness as a function of x 
can be adjusted. This choice reflects the fact that 
both the initial and final designs are characterized 
by symmetrical sections. The pressure distributions 
and airfoil sections for both the initial condition and 
the target are shown in Figure 1 1 . Figure 12 displays 
the solution after 19 design iterations. The target 
pressure distribution is almost obtained over most 
of the wing, with small discrepancies occurring close 
to the leading edge. There are no discernible differ- 
ences between the final airfoils and their targets. It 
is possible that there are not enough basis functions 
to allow exact recovery of the target pressure. How- 
ever, the discrepancies may also result simply from 
the fact that it would take a considerably greater 
number of design iterations to obtain better con- 
vergence to a more precise minimum. The design 
process was stopped when the rate of reduction in 
the cost function between design iterations slowed 
considerably. 

To explore the wing-body design capability of the 
general mesh formulation, a DC-9-30 planform is 
used as a testbed. The DC-9-30 is characterized 
by a straight tapered wing with an aspect ratio of 
8.7, \ chord sweep of 24.5°, washout of 4.5° and a 
taper ratio of 0.203. The fuselage has a diameter of 
1 1 .2 ft and a length of 107.6 ft. Figure 13 shows the 
surface mesh generated by WBGRID for the DC-9. 
Although the DC-9 cruises at a variety of different 
altitudes and Mach numbers its high speed cruise is 
M = 0.78. 

In our first attempt to redesign the wing con- 
tours for the DC-9-30, our starting point uses NACA 
64A410 airfoil defined at 9 stations across the span 
and scaled to 12.6 % thick at the root and 9.4 % 
thick at the tip. Even though we are using a rel- 
atively fine 161x41x41 mesh, a glance at Figure 13 
shows that the body mesh is still relatively coarse. 
But since body drag should stay essentially negligi- 
ble for inviscid calculations, the inverse of wing Lj D 
was used as the cost function. This amounts to wing 
optimization in the presence of the body. Both air- 
foil sections and pressure distributions for the ini- 
tial condition are shown in Figure 14. The NACA 
64A410 airfoil sections cause a strong shock across 
the entire upper surface of the wing since it was 
never intended for transonic use. A total of 1 17 de- 
sign variables are used to modify the wing. Nine de- 
sign variables adjust the twist at each wing defining 
station, while 54 design variables modify the leading 
edge shape and another 54 design variables alter the 
camber. The use of such a mix of variables demon- 
strates the flexibility of our choice of basis functions. 
The current choice fixes the planform and maximum 
thickness distribution but allows for a large range of 
wing shapes. Further, since the majority of the wing 
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is designed via variables that alter camber and twist, 
no prescribed pressure distributions were necessary 
to retain thickness. 

The wing was optimized in the constant Cl mode 
and after 16 design iterations it is shown in Figure 
15. The pressure distributions demonstrate that a 
significant reduction in the shock strength has been 
realized at the same time that the leading edge of the 
wing has been properly loaded. The wing drag coef- 
ficient has been reduced from 0.0248 to 0.0155 while 
the wing L/D has increased from 18.28 to 28.87. 
As in the case for the analytic mesh wing designs, 
a subsonic case was analyzed to determine the level 
of induced drag present in the configuration. These 
calculations revealed that the configuration contains 
149 counts of induced drag on the wing; thus the 
approximate wave drag coefficient for the optimized 
design is 0.0006. 

This calculation used 4.8 hours of Cray C-90 single 
processor time including all flow and adjoint calcu- 
lations. An estimate of the comparable calculation 
using finite differences for the gradient calculations 
is 160 hours. The calculation therefore shows a fac- 
tor of 33 savings in CPU time through the use of 
continuous sensitivities. 

A second attempt to redesign the DC-9 wing at 
a Mach number of 0.80 represents a more difficult 
challenge since it pushes the 24.5° swept wing well 
past its design conditions. Again 117 design vari- 
ables were specified such that twist, leading edge 
shape, and camber could be modified. The initial 
NACA 64A410 airfoil sections and pressure distri- 
butions are shown in Figure 16. Figure 17 shows 
the result after 15 design iterations where Cp w has 
been reduced from 0.0315 to 0.0167 and the wing lift 
to drag ratio has been increased from 14.44 to 26.86. 
The remaining component of the coefficient of drag 
due to shock waves for this design is about 0.0018. 
While this shows a dramatic performance improve- 
ment over the original wing it is still not quite as 
good as the wing designed at M = 0.78. This indi- 
cates that it may be difficult to obtain a very good 
design at this higher Mach number without allowing 
a change in thickness, a change in sweep, or a reduc- 
tion in the operating lift coefficient. An examination 
of the pressures reveals that this design is probably 
incurring the additional drag because of a stronger 
shock wave that traverses the span. 

Conclusions and Recommendations 

In the period since this approach to optimal shape 
design was first proposed by the author [18], the 
method has been verified by numerical implementa- 
tion for both potential flow and flows modeled by 
the Euler equations [19, 34, 23]. It has been demon- 
strated that it can be successfully used with a finite 
volume formulation to perform calculations with ar- 


bitrary numerically generated grids [34, 23]. Here 
results are presented for three-dimensional calcula- 
tions using both the analytic mapping and general 
finite volume implementations. The design of both 
wing and wing-body configurations indicates that 
this approach has matured to the extent that it can 
be a useful tool for the design of new aircraft. A 
factor of 33 savings in CPU time has been accom- 
plished through the use of the adjoint formulation. 
The clear limitation in these current results is the 
reliance on a single structured block for both the 
state and co-state fields. In the future our group in- 
tends not only to extend the method to treat both 
multiblock and unstructured meshes, but also to im- 
plement a Navier-Stokes version. 
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2a: Initial Wing 

() = 0.5001, C d = 0.0185, a = -0.958° 


2b: 8 Design Iterations 

C\ = 0.4929, C d = 0.0118, a = 0.172° 


Figure 2: Lifting Design Case, M = 0.82, Fixed Lift Mode. 
Drag Reduction 
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UPPER SURFACE PRESSURE 


LOWER SURFACE PRESSURE 


Figure 4: Lifting Design Case, M — 0.82, Fixed Lift Mode. 
Design after 8 cycles 
C L = 0.4929, C D = 0.01 18, « = 0.172° 

Drag Reduction 
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5a: span station z = 0.00 5b: span station z = 0.25 



5c: span station z = 0.50 5d: span station c = 0.75 


Figure 5 FL067 check on initial wing 

M = 0.82, C L = 0.4959, C D = 0.0171, a = -1.080° 
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6c: span station r = 0.50 6d: span station = 0.75 


Figure 6: FL067 check on redesigned wing. 

M = 0.82, C L = 0.4975, C D = 0.0107, a = 0.200° 
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7a: Initial Wing 

Ci = 0.5500, C d = 0.0243, a = -0.962° 


7b: 10 Design Iterations 

Ci = 0.5500, C d = 0.0144, « = 0.274° 


Figure 7: Lifting Design Case, M = 0.85, Fixed Lift Mode 
Drag Reduction 
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UPPER SURFACE PRESSURE 


UPPER SURFACE PRESSURE 


8a: Initial Wing 

Lifting Design Case, M = 0.85, Fixed Lift Mode 
Cl = 0.5500, C D = 0.0-243, a = -0.962° 
Drag Reduction 


8b: 10 Design Iterations 
Lifting Design Case, M = 0.85, Fixed Lift Mode 
C L = 0.5500, C D = 0.0144, a = 0.274° 
Drag Reduction 


Figure 8 Lifting Design Case, M - 0.85, Fixed Lift Mode. 
Drag Reduction 




9a: span station z = 0.00 



9b: span station z = 0.312 



9c: span station z = 0.625 9d: span station ^ = 0.937 


Figure 9: FL067 check on initial wing 

M = 0.85, C L = 0.5506, C D = 0.0236, a = -1.260° 


23 


ORKW4AL PAGE IS 

OF POOP QUALITY 




10c span station z - 0.625 10d: span station z - 0.937 

Figure 10: FL067 check on redesigned wing. 

M = 0 . 85 , C L = 0 . 5500 , C n = 0 . 0119 , a = 0 . 210 °. 
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11c: span station 2 = 0.604 lid: span station z = 0.896 


Figure 11: Initial condition for ONERA-M6 design. 

M = 0.84, C L = 0.1723, C D = 0.0122, a = 3.060° 

— , x Initial Wing: NACA 0012. 

,4- Target C p ONERA-M6 

150 Design Variables. 161x41x41 Mesh, Euler Wing/Body Solution 
Inverse Design 
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12a: span station z — 0.021 


12b: span station 2 = 0.312 
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12c span station z = 0.604 


1 2d : span station 2 = 0.896 


Figure 12: Solution after 19 iterations for ONERA-M6 design 
M = 0.84, C L = 0.1651, C D = 0.0078, a = 3.060°. 

— r x Initial Wing: NACA 0012. 

, + Target C p : ONERA-M6 

150 Design Variables, 161x41x41 Mesh, Euler Wing/Body Solution 
Inverse Design 
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Figure 13: DC-9-30 Planform and Wing-Body mesh 
161x41x41 C-H Topology 
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14a span station 2 = 5.85 ft. 14b: span station 2 = 18.02 ft 
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14c: span station z = 30.19 ft. 14d: span station r = 42.36 ft. 


Figure 14: Initial condition for DC-9-30 design, M = 0.78 
C L = 0.500, C Lm = 0.4535, C Dvt = 0.0248, c* = -3.769° 

117 Design Variables, 161x41x41 Mesh, Euler Wing/Body Solution 
Drag Minimization. 
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15a: span station z — 5.85 ft 



15b: span station z = 18.02 ft. 



15c: span station z = 30.19 ft. 15d: span station z = 42.36 ft 


Figure 15: DC-9-30 design after 16 iterations, M — 0.78 
C L = 0.500, C L „ = 0.4478, C Dv = 0.0155, a = -2.1776° 

117 Design Variables, 161x41x41 Mesh, Euler Wing/Body Solution 
Drag Minimization. 
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16a: span station z = 5.85 ft. 16b: span station z = 18.02 ft 



16c: span station z = 30.19 ft. 16d: span station z = 42.36 ft 


Figure 16: Initial condition for DC-9-30 design, M ~ 0.80 
C L = 0.500, C Lut — 0.4556, C Dvi = 0.0315, a = -3.960° 

117 Design Variables, 161x41x41 Mesh, Euler Wing/Body Solution 
Drag Minimization. 
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17a: span station z — 5.85 ft. 


17b: span station c = 18.02 ft 



17c: span station c = 30.19 ft 17d: span station z = 42.36 ft 


Figure 17: DC-9-30 design after 15 iterations M — 0.80 
C L = 0.500, C Lyt> = 0.4475, C Dw = 0.0167, a = -1.8374° 

117 Design Variables, 161x41x41 Mesh, Euler Wing/Body Solution 
Drag Minimization. 
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